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ABSTRACT 


Pseudo-bending  is  an  algorithm  for  calculating  seismic  travel  time  through  complex  3D  velocity  models.  The 
algorithm  was  originally  proposed  by  Um  and  Thurber  (1987)  and  later  extended  by  Zhao  et  al.  (1992)  to  account 
for  first  order  velocity  discontinuities.  We  have  modified  Zhao’s  method  of  handling  discontinuities  by 
implementing  a  two-dimensional  (2D)  minimization  algorithm  that  searches  for  the  point  on  the  velocity 
discontinuity  surface  where  Snell’s  Law  is  satisfied.  Further,  our  implementation  reduces  the  likelihood  that  the 
pseudo-bending  algorithm  will  return  a  local  minimum  by  starting  the  ray  calculation  from  several  different  starting 
rays.  Specifically,  interfaces  are  defined  that  include  first  order  discontinuities  plus  additional  interfaces  at  levels  of 
the  model  where  local  minima  might  be  generated.  Rays  are  computed  that  are  constrained  to  bottom  in  each  layer 
between  these  interfaces.  The  computed  rays  might  be  reflected  off  the  top  of  the  layer,  turn  within  the  layer,  or 
diffract  along  the  interface  at  the  bottom  of  the  layer.  The  computed  ray  that  is  seismologically  valid  and  that  has  the 
shortest  travel  time  is  retained. 

The  modifications  we  have  made  to  the  algorithm  have  made  it  more  accurate  and  robust  but  have  also  made  it  more 
computationally  expensive.  To  mitigate  this  impact,  we  have  implemented  our  software  in  a  distributed  parallel 
computing  environment,  which  makes  possible  the  calculation  of  many  rays  simultaneously. 
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OBJECTIVE 


Our  goal  is  to  develop  a  practical,  accurate,  infinite  frequency  travel-time  calculator  through  three-dimensional  (3D) 
models  of  the  distribution  of  seismic  velocity  in  the  Earth. 

RESEARCH  ACCOMPLISHED 
Background 

Three-dimensional  models  of  the  velocity  distribution  within  the  Earth  are  becoming  increasingly  available  in  the 
seismological  research  community  and  have  the  potential  to  significantly  improve  our  ability  to  accurately  and 
precisely  locate  seismic  events  around  the  world  (e.g,.  Flanagan  et  ah,  2007).  These  3D  models  are  generally  derived 
directly  from  observed  travel  times,  using  tomographic  inversion  techniques.  Development  of  these  models  and  use 
of  the  models  in  seismic-event  location  calculations,  requires  the  ability  to  accurately  compute  predicted  source-to- 
receiver  travel  times  though  the  proposed  3D  velocity  structures.  This  is  a  complicated  and  computationally 
expensive  task,  and  the  many  algorithms  that  have  been  proposed  to  accomplish  it  fall  into  three  broad  categories. 
Eikonal  solvers  (e.g.,  Vidale,  1988,  1990;  Podvin  and  Lecomte,  1991;  Rawlinson  and  Sambridge,  2004a,  2004b;  and 
deKool  et  al.,  2006)  employ  finite-difference  techniques  to  propagate  wavefronts  through  3D  velocity  distributions 
defined  on  a  3D  grid  of  points,  thereby  computing  predicted  travel  times  between  a  single  point  and  an  entire  grid  of 
points  distributed  within  the  medium.  Ray  shooters  (e.g.,  Menke)  systematically  perturb  an  initial  estimate  of  the  ray 
take-off  angle  from  a  seismic  source  until  the  ray  hits  the  seismic  receiver,  within  some  prescribed  tolerance.  Ray 
benders  (e.g.,  Um  and  Thurber,  1987)  start  with  an  initial  estimate  of  the  ray  geometry  connecting  the  source  and 
receiver  and  perturb  the  geometry  of  the  ray  until  Fermat’s  principle  of  stationary  time  is  satisfied  all  along  the  ray 
path.  Though  very  different  in  approach  to  the  problem,  all  of  these  techniques  assume  that  the  wave  propagates 
from  source  to  receiver  with  infinite  frequency. 

After  reviewing  several  available  algorithms,  we  chose  to  implement  a  version  of  the  ray-bending  algorithm 
described  by  Um  and  Thurber  (1987),  which  was  extended  to  include  the  effects  of  velocity  discontinuities  by  Zhao 
et  al.  (1992)  and  Koketsu  and  Sekine  (1998).  We  considered  this  algorithm  superior  to  finite-difference  eiknonal 
solvers  for  our  intended  applications  because  it  provides  direct,  point-to-point  travel  times  without  interpolation  off 
a  grid  and  because  it  does  not  require  calculation  of  full  volumes  of  travel  times  for  a  small  number  of 
source-receiver  pairs.  We  also  consider  it  to  be  superior  to  ray  shooters  in  that  it  will  always  find  a  solution,  even  in 
the  presence  of  low- velocity  zones,  which  can  cause  ray  shooters  to  fail. 

The  algorithm  is  based  on  the  fact  that,  along  a  valid  ray 
path,  the  ray  curvature  is  proportional  to  the  component 
of  the  local  velocity  gradient  normal  to  the  ray  path 
(Figure  1).  The  algorithm  steps  along  an  initial  estimate 
of  the  path,  perturbing  the  center  point  (Xk  in  Figure  1)  of 
three  adjacent  points  on  the  path  so  that  it  satisfies  this 
condition  (Xk'  in  Figure  1).  Zhao  et  al.  (1992)  and 
Koketsu  and  Sekine  (1998)  extended  the  algorithm  to 
include  the  effects  of  planar  and  spherical  velocity 
discontinuities,  respectively.  In  this  paper,  we  further 
extend  the  algorithm  to  allow  discontinuities  at  arbitrary 
orientations.  These  modifications,  described  in  detail  in 
the  next  section,  preclude  the  inclusion  of  the 
“enhancement  factor”  described  by  Um  and  Thurber 
(1987),  and  we  have  not  included  it  in  our  implementation 

A  key  feature  of  any  system  for  computing  travel  time  through  3D  Earth  models  is  the  subsystem  used  to  represent, 
and  interpolate  velocity  values  from,  the  3D  model  itself.  The  system  we  use  is  described  in  detail  in  a  companion 
paper  (Ballard  et  al.,  2008;  these  Proceedings)  and  is  not  reviewed  here. 


Figure  1.  Pseudo-bending  algorithm  (from  Um 
and  Thurber,  1987). 
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Velocity  Discontinuities 

Figure  2  depicts  two  media,  Mi  and  M2,  with  velocity  distributions  V/(x,  y,  z)  and  V2(x,  y,  z),  which  may  vary  in 
three  dimensions  but  which  are  first-order  continuous  internally.  Surface  S  represents  the  smooth  boundary  between 
the  two  media  such  that  the  velocity  may  be  discontinuous  across  S.  The  curve  in  Figure  2  represents  the 
intersection  of  S  with  the  plane  of  the  figure,  but  S  may  dip  into  and/or  out  of  the  plane  of  the  figure.  Unit  vector  n 
is  normal  to  S  at  X  and  does  not  necessarily  lie  in  the  plane  of  the  figure.  Pi  and  P2  are  points  in  Mi  and  M2, 
respectively.;  xx  and  x2  are  unit  vectors  that  point  from  X  toward  Pi  and  P2  respectively. 


In  order  for  points  Pi,  X  and  P2  to  lie  on  a  common  ray  path,  Snell’s  Law  must  be  honored  at  X,  which  requires  that 
the  following  3  conditions  be  satisfied: 


sin  6^ 

=_U 


(i) 


C2  =  (Xj  x  x2 )  •  n  =  0 


(2) 


and 

C3  =  (xj  x  n) •  (x2  x  n)  <  0  .  (3) 

In  Equation  l,  Vi,  z=l  ,2  is  the  arithmetic  average  of  the  velocity  at  Pi  and  the  velocity  at  X  on  the  same  side  of  S  as 
Pi.  Equation  2  requires  that,  xx ,  x2 ,  and  n  be  coplanar,  implying  that  n  lies  in  the  plane  of  Figure  2  when  Snell’s 
Law  is  satisfied.  Equation  3  prohibits  xx  and  x2  from  residing  on  the  same  side  of  the  line  containing  n  (i.e.,  a 
nonphysical  ray  path).  Figure  3  (a,b,  and  c)  illustrates  example  distributions  of  Equations  1  through  3  displayed  on  a 
projection  of  S  onto  a  spherical  cap  that  passes  through  X. 

While  ray  tracing,  we  are  given  points  Pi  and  P2  on  opposite  sides  of  S  and  seek  point  X  on  S  that  satisfies  Snell’s 
Law,  as  expressed  by  the  three  conditions  above.  To  find  X,  we  implement  a  2D  search  algorithm  that  minimizes 

/•  =  e  -  |c,  ~ioc;|.  (4) 

Experience  indicates  that  the  factor  of  10  applied  to  C3  improves  convergence  behavior.  Figure  3d  illustrates  an 
example  distribution  of  F  on  a  projection  of  S  onto  a  spherical  cap  that  passes  through  X. 


Pi 


Figure  2.  Snell’s  Law  applied  at  point  X  on 
surface  S. 
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Phase  Discrimination 

A  problem  that  frequently  arises  when  using  the  bending  algorithm  to  predict  travel  times  through  realistic,  regional 
to  teleseismic  scale,  3D  models  is  that  there  exists  more  than  one  valid  ray  that  connects  a  given  source-receiver  pair 
(Figure  4).  Some  of  these  rays  represent  different  seismic  phases,  but  even  for  a  single  phase  there  can  be  multiple 
valid  rays  that  bottom  at  different  levels  in  the  Earth.  Each  of  these  represents  a  local  minimum  travel  time  but  only 
once  can  represent  the  overall  minimum  path.  To  ensure  that  we  get  the  desired  ray,  we  specify  a  number  of 
interfaces  within  the  Earth  and  require  the  bender  to  compute  rays  that  are  constrained  to  bottom  in  each  of  the 
layers  defined  by  those  interfaces.  We  specify  interfaces  at  all  the  velocity  discontinuities  with  which  a  ray  of  a 
given  phase  must  interact  (Moho,  410,  660,  etc.).  We  specify  additional  interfaces  at  levels  within  the  Earth  model 
that  could  potentially  generate  travel-time  triplications,  even  though  the  velocity  is  continuous  across  these 
interfaces.  The  ray  computed  for  each  layer  might  be  one  that  refracts  within  the  layer,  reflects  off  the  interface  at 
the  top  of  the  layer,  or  diffracts  along  the  interfaces  at  the  top  and/or  bottom  of  the  layer. 
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Figure  3.  Plots  of  (a)  Equation  1,  (b)  Equation  2,  (c)  Equation  3,  and  (d)  Equation  4  displayed  on  a 

projection  of  S  onto  a  spherical  cap  that  passes  through  X.  Points  PI  and  P2  lie  above  and  below 
the  plane  of  the  figures,  respectively.  In  (d),  portions  of  the  map  where  Equation  3  is  not  equal  to 
zero  are  blanked  out. 
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Not  all  of  these  rays  are  valid  and  hence  can  be  discarded.  In  particular,  rays  that  reflect  off  of,  or  diffract  along, 
defined  interfaces  across  which  the  velocity  is  continuous  are  not  valid.  If  more  than  one  valid  ray  is  calculated  from 
a  given  set  of  defined  interfaces,  the  ray  with  the  shortest  travel  time  is  retained. 

A  further  advantage  of  this  approach  is  that  travel  time  for  any  desired  phase  can  be  computed  by  specifying  the 
interfaces  in  the  Earth  with  which  the  ray  must  interact  (Figure  5).  For  example,  to  compute  predictions  for  phase 
Pn,  one  would  specify  only  the  Moho  and  other  interfaces  in  the  mantle  down  to,  but  not  including,  the  660-km 
discontinuity,  thereby  constraining  the  computed  ray  to  bottom  somewhere  in  the  mantle  above  the  660. 
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Figure  4.  Multiple  valid  seismic  rays  connect  a  given  source  and  receiver,  each  constrained  to  bottom  in  a 
specified  layer  of  the  Earth  model. 
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Figure  5.  Travel-time  curves  for  a  source-receiver  pair  in  a  3D  Earth  model  with  velocity  generally  less 
than  that  of  the  AK135  model. 


Algorithm  Verification 

We  have  verified  the  validity  of  the  bender  algorithm  in  two  ways.  First,  we  constructed  a  3D  version  of  the  radially 
symmetric  AK135  model  (Kennett  et  al.,  1995).  The  3D  version  is  deployed  on  a  uniform,  1°  triangular  tessellation 
similar  to  the  tessellations  we  use  for  more-complex  3D  models  (Ballard  et  al.,  2008;  these  Proceedings),  but  the 
model  is  based  on  a  spherical  Earth  with  the  same  velocity  discontinuities  and  radial  velocity  distribution  as  the 
AK135  model.  Rays  through  the  model  computed  with  Bender  are  shown  in  Figure  6,  and  travel  times  computed 
with  Bender  and  using  the  TaupToolkit  software  (Crotwell  et  al.,  1999)  are  compared  in  Figure  7.  The  difference  at 
all  regional  and  teleseismic  distances  is  at  most  a  few  hundredths  of  seconds,  indicating  that  Bender  is  very  accurate, 
even  for  long  path  lengths. 


Figure  6.  Rays  computed  with  Bender  through  a  3D  version  of  the  AK135  Earth  model. 
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To  verify  the  validity  of  Bender  travel  times 
through  Earth  models  with  significant  3D  velocity 
variations,  we  built  a  3D  model  based  on  the  a 
priori  GNEMRD  Unified  Model  (Begnaud  et  ah, 

2004;  Pasyanos  et  ah,  2004)  and  computed  travel 
times  through  it  using  our  Bender  code  and  the 
FMM  code  (De  Kool  et  ah,  2006).  Topography  and 
bathymetry  for  the  Unified  Model  were  derived 
from  ETOP05,  while  shallow  sedimentary  layers 
were  taken  from  the  1  -degree-resolution  sediment 
model  of  Laske  and  Masters  (1997).  The  structure 
of  crust  and  uppermost  mantle  is  drawn  from  a 
variety  of  focused  regional  studies.  The 
upper-mantle  model  (below  100  km)  is  taken  from 
the  a  priori  3SMAC  model  (Nataf  and  Ricard, 

1996). 

As  with  all  finite-difference  eikonal  solvers,  the 
FMM  algorithm  works  by  following  a  wavefront  as 
it  moves  across  a  volume  of  grid  points,  updating 
the  travel  times  in  the  grid  according  to  the  eikonal 
differential  equation,  using  a  second-order  finite-difference  scheme.  We  chose  to  use  FMM  for  our  comparison 
because  it  has  a  number  of  important  features  that  make  it  well  suited  for  making  calculations  in  the 
regional-to-teleseismic  distance  ranges  that  we  are  interested  in.  First,  the  propagation  grid  is  in  spherical 
coordinates  rather  than  Cartesian  coordinates,  so  the  sphericity  of  the  Earth  can  be  more  readily  represented.  Second, 
interfaces  are  explicitly  represented  in  the  FMM  grid,  regardless  of  the  propagation  grid  spacing.  Thus,  a  finer 
propagation  grid  is  not  necessary  to  capture  discontinuities  with  topography,  like  the  Moho.  Third,  FMM  allows 
calculation  of  secondary  phases,  so  it  can  be  used  for  non-first- arrival  phases.  Finally,  finer  grid  spacing  can  be  used 
near  the  source  to  better  capture  the  highly  curved  wavefront  and  hence  improve  accuracy  without  drastically 
increasing  the  overall  number  of  nodes. 

For  purposes  of  comparing  travel  times  computed  using  Bender  with  those  computed  using  FMM,  we  chose  station 
NIF  located  in  Nilore,  Pakistan,  at  latitude  33.6506  N  and  longitude  72.2686  E.  Even  though  the  station  is  0.629  km 
above  sea  level,  we  set  the  elevation  to  sea  level  for  these  calculations  in  order  to  avoid  having  to  compute  elevation 
corrections  for  the  FMM  results.  Next,  we  specified  a  40°  x  40°  grid  of  sources  extending  from  53°  to  93°  E 
longitude  and  14°  to  54°  N  latitude,  with  0.25°  grid  spacing  in  both  directions.  All  sources  were  positioned  at  10  km 
below  sea  level.  We  computed  travel  times  from  these  sources  to  the  position  of  station  NIF  using  both  Bender  and 
FMM.  With  Bender,  it  was  possible  to  compute  point-to-point  travel  times  for  these  source-receiver  pairs  directly, 
but  with  FMM  we  had  to  first  construct  a  3D  grid  that  extended  from  sea  level  to  a  depth  of  1000  km  in  order  to 
ensure  that  the  depth  of  the  grid  exceeded  the  deepest  bottoming  depth  of  any  rays  we  were  interested  in.  In 
addition,  the  FMM  grid  had  to  be  denser  than  0.25°  x  0.25°  in  order  to  achieve  sufficient  accuracy.  To  test  the 
accuracy  of  the  FMM  results,  we  constructed  two  finite-difference  grids,  one  0.125°  x  0.125°  x  10  km  and  the  other 
0.0625°  x  0.0625°  x  5  km.  Due  to  computer-memory  limitations,  we  had  to  break  the  3D  grid  into  four  quadrants, 
run  FMM  on  each  quadrant  separately,  recombine  the  results  into  a  single  grid,  and  finally  subsample  the  FMM 
grids  at  the  points  of  the  2D  Bender  grid.  The  difference  in  travel  time  between  FMM  and  Bender,  for  the  two 
different  FMM  grid  resolutions,  is  illustrated  in  Figure  8.  The  difference  in  travel  time  decreased  as  a  function  of 
FMM  grid  resolution,  and  FMM  travel  time  was  virtually  always  greater  than  Bender  travel  time.  From  these 
observations  we  conclude  that  the  Bender  travel  times  are  likely  a  better  reflection  of  the  actual  travel  time  through 
the  Earth  model  used  in  the  calculations.  It  is  entirely  possible  that  FMM  travel  times  would  have  been  reduced,  and 
approached  the  Bender  travel  times,  had  we  been  able  to  improve  the  resolution  of  the  FMM  grid  even  further.  We 
were  prevented  from  doing  that  by  computer-memory  limitations 


Figure  7.  Difference  in  travel  time  computed  by  Bender 
and  Taup  Toolkit  through  the  AK135 
velocity  model. 
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FMM  Coarse  Resolution 


FMM  Fine  Resolution 


Figure  8.  FMM  travel  time  minus  Bender  travel  time  for  a  40°  x  40°  grid  surrounding  station  NIL.  On  the 
left,  the  FMM  travel  times  were  computed  on  a  0.125°  x  0.125°  x  10  km  grid;  on  the  right,  they 
were  computed  on  a  0.0625°  x  0.0625°  x  5  km  grid.  Black  contour  lines  indicate  positions  of 
sources  where  rays  bottomed  at  depths  of  100  km,  440  km,  and  660  km,  as  determined  by 
Bender. 


Distributed  Computing 

The  enhancements  we  have  made  to  the  pseudo-bending  algorithm  of  Um  and  Thurber  (1987),  as  modified  by  Zhao 
et  al.  (1992),  have  made  it  more  robust  and  accurate  but  have  also  made  it  more  computationally  expensive.  Our 
implementation  typically  requires  approximately  0.5  seconds  to  compute  one  travel  time  on  a  single  computer 
processor.  To  mitigate  these  negative  performance  impacts,  we  have  implemented  the  software,  which  is  written 
entirely  in  the  Java  programming  language,  in  a  distributed  computing  environment,  using  the  Java  Parallel 
Processing  Framework  (JPPF;  www.ippf.org).  This  technology  makes  it  fairly  straightforward  to  run  many 
travel-time  calculations  simultaneously  on  all  available  PC,  Mac,  Unix,  and  Linux  computers.  We  are  currently 
running  our  software  on  a  cluster  of  nine  16-processor  PC  computers  running  the  64-bit  Windows  Server  2003 
operating  system.  Combined  with  other  PCs,  Linux  boxes,  and  Sun  workstations,  our  system  has  almost  200 
processors  available  for  calculating  travel  times. 


CONCLUSIONS 

We  have  implemented  the  pseudo-bending  algorithm  of  Um  and  Thurber  (1987),  as  modified  by  Zhao  et  al.  (1992), 
for  computing  predicted  travel  times  through  3D  velocity  models.  We  have  modified  Zhao’s  method  of  handling 
discontinuities  by  implementing  a  2D  minimization  algorithm  that  searches  for  the  point  on  the  velocity 
discontinuity  surface  where  Snell’s  Law  is  satisfied.  Our  implementation  also  significantly  decreases  the  likelihood 
that  the  algorithm  will  return  a  ray  that  represents  a  local  minimum  by  computing  rays  that  bottom  in  many  different 
specified  layers  in  the  model  and  returning  the  ray  that  yields  the  smallest  travel  time.  We  have  validated  our 
algorithm  by  comparing  it  with  the  TaupToolkit  software  (Crotwell  et  al.,  1999)  for  radially  symmetric  ID  Earth 
models.  We  have  also  compared  our  results  with  travel  times  computed  with  the  FMM  finite-difference  calculator 
(De  Kool  et  al.,  2006)  and  demonstrated  that  the  Bender  yields  superior  results  with  substantially  fewer  computer 
resources. 

We  are  currently  developing  improved  velocity  models  of  the  Earth  using  a  tomographic  inversion  system  that  uses 
our  Bender  software  as  the  forward  travel-time  calculator.  Preliminary  results  for  a  study  in  southcentral  Asia  are 
reported  in  Young  et  al.  (2008,  these  Proceedings). 
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